#' ---
#' title: "Robustness Checks 1 and 2"
#' author: "Gento Kato & Fan Lu"
#' date: "June 8, 2023"
#' ---
#' 
#' 
#' # Preparation 
#' 

## Clean Up Space
rm(list=ls())

## Set Working Directory (Automatically) ##
setwd(dirname(rstudioapi::getActiveDocumentContext()$path)); 

require(interflex)
require(estimatr)
require(texreg)
require(ggplot2)
require(ggrepel)
require(haven)

survey22 <- readRDS("data_survey22_v7.rds")
icpp09 <- readRDS("data_icpp09_v7.rds")
icpp13 <- readRDS("data_icpp13_v7.rds")

cfmap <- list("(Intercept)"="(Intercept)",
              "edu2"="University Education",
              "cohortCohort II (18+ in 1976-1989)"="Cohort II (18+ in 1975-89)",
              "cohortCohort III (18+ in 1990-99)"="Cohort III (18+ in 1990-99)",
              "cohortCohort IV (18+ in 2000-09)"="Cohort IV (18+ in 2000-09)",
              "cohortCohort V (18+ in 2010-)"="Cohort V (18+ in 2010-)",
              "edu2:cohortCohort II (18+ in 1976-1989)"="University * Cohort II",
              "edu2:cohortCohort III (18+ in 1990-99)"="University * Cohort III",
              "edu2:cohortCohort IV (18+ in 2000-09)"="University * Cohort IV",
              "edu2:cohortCohort V (18+ in 2010-)"="University * Cohort V",
              "male"="Gender (Male)", 
              "cohortCohort II (18+ in 1976-1989):male"="Male * Cohort II",
              "cohortCohort III (18+ in 1990-99):male"="Male * Cohort III",
              "cohortCohort IV (18+ in 2000-09):male"="Male * Cohort IV",
              "cohortCohort V (18+ in 2010-):male"="Male * Cohort V",
              "incomecatMiddle" = "Income (Middle)",
              "cohortCohort II (18+ in 1976-1989):incomecatMiddle"="Income (Middle) * Cohort II",
              "cohortCohort III (18+ in 1990-99):incomecatMiddle"="Income (Middle) * Cohort III",
              "cohortCohort IV (18+ in 2000-09):incomecatMiddle"="Income (Middle) * Cohort IV",
              "cohortCohort V (18+ in 2010-):incomecatMiddle"="Income (Middle) * Cohort V",
              "incomecatHigh" = "Income (High)",
              "cohortCohort II (18+ in 1976-1989):incomecatHigh"="Income (High) * Cohort II",
              "cohortCohort III (18+ in 1990-99):incomecatHigh"="Income (High) * Cohort III",
              "cohortCohort IV (18+ in 2000-09):incomecatHigh"="Income (High) * Cohort IV",
              "cohortCohort V (18+ in 2010-):incomecatHigh"="Income (High) * Cohort V",
              "incomecatMissing" = "Income (Missing)",
              "cohortCohort II (18+ in 1976-1989):incomecatMissing"="Income (Missing) * Cohort II",
              "cohortCohort III (18+ in 1990-99):incomecatMissing"="Income (Missing) * Cohort III",
              "cohortCohort IV (18+ in 2000-09):incomecatMissing"="Income (Missing) * Cohort IV",
              "cohortCohort V (18+ in 2010-):incomecatMissing"="Income (Missing) * Cohort V",
              "workstatStudent/Housemaker/Part-Time" = 
                "Student/Housemaker/Part-Time",
              "cohortCohort II (18+ in 1976-1989):workstatStudent/Housemaker/Part-Time"=
                "Student/Housemaker/Part-Time * Cohort II",
              "cohortCohort III (18+ in 1990-99):workstatStudent/Housemaker/Part-Time"=
                "Student/Housemaker/Part-Time * Cohort III",
              "cohortCohort IV (18+ in 2000-09):workstatStudent/Housemaker/Part-Time"=
                "Student/Housemaker/Part-Time * Cohort IV",
              "cohortCohort V (18+ in 2010-):workstatStudent/Housemaker/Part-Time"=
                "Student/Housemaker/Part-Time * Cohort V",
              "workstatStudent/Part-Time" = 
                "Student/Part-Time",
              "cohortCohort II (18+ in 1976-1989):workstatStudent/Part-Time"=
                "Student/Part-Time * Cohort II",
              "cohortCohort III (18+ in 1990-99):workstatStudent/Part-Time"=
                "Student/Part-Time * Cohort III",
              "cohortCohort IV (18+ in 2000-09):workstatStudent/Part-Time"=
                "Student/Part-Time * Cohort IV",
              "cohortCohort V (18+ in 2010-):workstatStudent/Part-Time"=
                "Student/Part-Time * Cohort V",
              "workstatSelf-Employed/Full-Time/Managerial" = 
                "Self-Employed/Full-Time",
              "cohortCohort II (18+ in 1976-1989):workstatSelf-Employed/Full-Time/Managerial"=
                "Self-Employed/Full-Time * Cohort II",
              "cohortCohort III (18+ in 1990-99):workstatSelf-Employed/Full-Time/Managerial"=
                "Self-Employed/Full-Time * Cohort III",
              "cohortCohort IV (18+ in 2000-09):workstatSelf-Employed/Full-Time/Managerial"=
                "Self-Employed/Full-Time * Cohort IV",
              "cohortCohort V (18+ in 2010-):workstatSelf-Employed/Full-Time/Managerial"=
                "Self-Employed/Full-Time * Cohort V",
              "married"="Currently Married",
              "cohortCohort II (18+ in 1976-1989):married"="Married * Cohort II",
              "cohortCohort III (18+ in 1990-99):married"="Married * Cohort III",
              "cohortCohort IV (18+ in 2000-09):married"="Married * Cohort IV",
              "cohortCohort V (18+ in 2010-):married"="Married * Cohort V",
              "urban"="Urbanness (Current Residence)",
              "cohortCohort II (18+ in 1976-1989):urban"="Urbanness * Cohort II",
              "cohortCohort III (18+ in 1990-99):urban"="Urbanness * Cohort III",
              "cohortCohort IV (18+ in 2000-09):urban"="Urbanness * Cohort IV",
              "cohortCohort V (18+ in 2010-):urban"="Urbanness * Cohort V")

#'
#' # Analysis
#'

#+
###################################################
## Analysis of 2009 Survey (Figure A1, Table A2) ##
###################################################

## Models with Categorical Interaction
mA_icpp09 <- lm_robust(immigincrease_alt~edu2*cohort+
                         male*cohort+
                         incomecat*cohort+
                         workstat*cohort+
                         married*cohort+
                         urban*cohort, 
                       data=icpp09, se_type="stata", clusters=as.factor(area))
mB_icpp09 <- lm_robust(foreignrights~edu2*cohort+
                         male*cohort+
                         incomecat*cohort+
                         workstat*cohort+
                         married*cohort+
                         urban*cohort, 
                       data=icpp09, se_type="stata", clusters=as.factor(area))
screenreg(list(mA_icpp09,mB_icpp09), 
          digits=3, include.ci=FALSE, single.row=TRUE,
          custom.model.names = c("Admit More Foreigners","Integrate Foreign Residents"),
          custom.coef.map = cfmap)
texreg(list(mA_icpp09,mB_icpp09), 
       digits=3, include.ci=FALSE, single.row=TRUE,
       custom.model.names = c("Admit More Foreigners","Integrate Foreign Residents"),
       custom.coef.map = cfmap, 
       custom.note = "%stars. Robust standard errors in parentheses.",
       booktabs = TRUE, dcolumn = TRUE, use.packages = FALSE, table=FALSE,
       file = "mtab_rc1_icpp09.tex")

mA_edu2_icpp09 <- 
  as.data.frame(  rbind(
    summary(lm_robust(immigincrease_alt~edu2*relevel(cohort,levels(cohort)[1])+
                        male*relevel(cohort,levels(cohort)[1])+
                        incomecat*relevel(cohort,levels(cohort)[1])+
                        workstat*relevel(cohort,levels(cohort)[1])+
                        married*relevel(cohort,levels(cohort)[1])+
                        urban*relevel(cohort,levels(cohort)[1]), 
                      data=icpp09, 
                      se_type="stata", clusters=as.factor(area)))$coefficients[2,c(1,4,5,6)],
    summary(lm_robust(immigincrease_alt~edu2*relevel(cohort,levels(cohort)[2])+
                        male*relevel(cohort,levels(cohort)[2])+
                        incomecat*relevel(cohort,levels(cohort)[2])+
                        workstat*relevel(cohort,levels(cohort)[2])+
                        married*relevel(cohort,levels(cohort)[2])+
                        urban*relevel(cohort,levels(cohort)[2]), 
                      data=icpp09, 
                      se_type="stata", clusters=as.factor(area)))$coefficients[2,c(1,4,5,6)],
    summary(lm_robust(immigincrease_alt~edu2*relevel(cohort,levels(cohort)[3])+
                        male*relevel(cohort,levels(cohort)[3])+
                        incomecat*relevel(cohort,levels(cohort)[3])+
                        workstat*relevel(cohort,levels(cohort)[3])+
                        married*relevel(cohort,levels(cohort)[3])+
                        urban*relevel(cohort,levels(cohort)[3]), 
                      data=icpp09, 
                      se_type="stata", clusters=as.factor(area)))$coefficients[2,c(1,4,5,6)],
    summary(lm_robust(immigincrease_alt~edu2*relevel(cohort,levels(cohort)[4])+
                        male*relevel(cohort,levels(cohort)[4])+
                        incomecat*relevel(cohort,levels(cohort)[4])+
                        workstat*relevel(cohort,levels(cohort)[4])+
                        married*relevel(cohort,levels(cohort)[4])+
                        urban*relevel(cohort,levels(cohort)[4]), 
                      data=icpp09, 
                      se_type="stata", clusters=as.factor(area)))$coefficients[2,c(1,4,5,6)]
  )  )
colnames(mA_edu2_icpp09) <- c("est","p","lci","uci")
mA_edu2_icpp09$dv <- "Admit More Foreigners"
# unique(sort(icpp09$univyr))
# median(1948:1974); median(1975:1989); median(1990:1999); median(2000:2008) 
mA_edu2_icpp09$x <- c(1961,1982,1994.5,2004)

mB_edu2_icpp09 <- 
  as.data.frame(  rbind(
    summary(lm_robust(foreignrights~edu2*relevel(cohort,levels(cohort)[1])+
                        male*relevel(cohort,levels(cohort)[1])+
                        incomecat*relevel(cohort,levels(cohort)[1])+
                        workstat*relevel(cohort,levels(cohort)[1])+
                        married*relevel(cohort,levels(cohort)[1])+
                        urban*relevel(cohort,levels(cohort)[1]), 
                      data=icpp09, 
                      se_type="stata", clusters=as.factor(area)))$coefficients[2,c(1,4,5,6)],
    summary(lm_robust(foreignrights~edu2*relevel(cohort,levels(cohort)[2])+
                        male*relevel(cohort,levels(cohort)[2])+
                        incomecat*relevel(cohort,levels(cohort)[2])+
                        workstat*relevel(cohort,levels(cohort)[2])+
                        married*relevel(cohort,levels(cohort)[2])+
                        urban*relevel(cohort,levels(cohort)[2]), 
                      data=icpp09, 
                      se_type="stata", clusters=as.factor(area)))$coefficients[2,c(1,4,5,6)],
    summary(lm_robust(foreignrights~edu2*relevel(cohort,levels(cohort)[3])+
                        male*relevel(cohort,levels(cohort)[3])+
                        incomecat*relevel(cohort,levels(cohort)[3])+
                        workstat*relevel(cohort,levels(cohort)[3])+
                        married*relevel(cohort,levels(cohort)[3])+
                        urban*relevel(cohort,levels(cohort)[3]), 
                      data=icpp09, 
                      se_type="stata", clusters=as.factor(area)))$coefficients[2,c(1,4,5,6)],
    summary(lm_robust(foreignrights~edu2*relevel(cohort,levels(cohort)[4])+
                        male*relevel(cohort,levels(cohort)[4])+
                        incomecat*relevel(cohort,levels(cohort)[4])+
                        workstat*relevel(cohort,levels(cohort)[4])+
                        married*relevel(cohort,levels(cohort)[4])+
                        urban*relevel(cohort,levels(cohort)[4]), 
                      data=icpp09, 
                      se_type="stata", clusters=as.factor(area)))$coefficients[2,c(1,4,5,6)]
  )  )
colnames(mB_edu2_icpp09) <- c("est","p","lci","uci")
mB_edu2_icpp09$dv <- "Integrate Foreign Residents"
mB_edu2_icpp09$x <- c(1961,1982,1994.5,2004)

m_edu2_icpp09 <- rbind(mA_edu2_icpp09,mB_edu2_icpp09)
m_edu2_icpp09$dv <- factor(m_edu2_icpp09$dv,levels=unique(m_edu2_icpp09$dv))
m_edu2_icpp09$lb <- factor(c("I","II","III","IV"), levels=c("I","II","III","IV"))
m_edu2_icpp09

## Interflex (without background variables)
set.seed(12345)
mxA_icpp09 <- interflex(Y = "immigincrease_alt", D = "edu2", X = "univyr",
                        Z = c("male", "incomecat", "workstat", "married", "urban"),
                        cl="area",
                        data = icpp09, na.rm=TRUE, estimator = "kernel",
                        full.moderate = TRUE, bw=5, #CI=FALSE,
                        nboots = 1000, parallel = TRUE, cores = 3, theme.bw=TRUE)
# mxA_icpp09$figure
set.seed(12345)
mxB_icpp09 <- interflex(Y = "foreignrights", D = "edu2", X = "univyr",
                        Z = c("male", "incomecat", "workstat", "married", "urban"),
                        cl="area",
                        data = icpp09, na.rm=TRUE, estimator = "kernel",
                        full.moderate = TRUE, bw=5, #CI=FALSE,
                        nboots = 1000, parallel = TRUE, cores = 3, theme.bw=TRUE)
# mxB_icpp09$figure

# save(mxA_icpp09, mxB_icpp09, file="mx_icpp09_alt.RData")
# load("mx_icpp09_alt.RData")

tmpA <- as.data.frame(mxA_icpp09$est.kernel$`1`)
tmpA$dv <- "Admit More Foreigners"
tmpB <- as.data.frame(mxB_icpp09$est.kernel$`1`)
tmpB$dv <- "Integrate Foreign Residents"
out <- rbind(tmpA,tmpB)
out$dv <- factor(out$dv,levels=unique(out$dv))

tmpA <- data.frame(x=mxA_icpp09$hist.out$mids,
                   y=c(mxA_icpp09$count.tr$`0`,mxA_icpp09$count.tr$`1`),
                   tr=rep(c("0","1"), each=length(mxA_icpp09$hist.out$mids)))
tmpA$dv <- "Admit More Foreigners"
tmpB <- data.frame(x=mxB_icpp09$hist.out$mids,
                   y=c(mxB_icpp09$count.tr$`0`,mxB_icpp09$count.tr$`1`),
                   tr=rep(c("0","1"), each=length(mxB_icpp09$hist.out$mids)))
tmpB$dv <- "Integrate Foreign Residents"
outh <- rbind(tmpA,tmpB)
outh$dv <- factor(outh$dv,levels=unique(outh$dv))
outh$tr <- factor(outh$tr, levels=c("1","0"))

#+ fig.width=7, fig.height=4
## Figure
movement0 <- 0.3
ggplot(out) + 
  geom_vline(aes(xintercept=1975), linetype=1, color="gray80") + 
  geom_vline(aes(xintercept=1990), linetype=1, color="gray80") + 
  geom_vline(aes(xintercept=2000), linetype=1, color="gray80") + 
  geom_vline(aes(xintercept=2010), linetype=1, color="gray80") + 
  geom_col(data=outh, aes(x=x,y=y*0.0015,fill=tr), width=0.8) + 
  geom_hline(aes(yintercept=0+movement0), linetype=2, color="gray50") + 
  geom_ribbon(aes(x=X,ymin=`lower CI(95%)`+movement0,
                  ymax=`upper CI(95%)`+movement0),alpha=0.3) +
  geom_line(aes(x=X,y=TE+movement0)) +
  # geom_line(aes(x=X, y=`lower CI(95%)`+movement0),alpha=0.7, linetype=2) +
  # geom_line(aes(x=X, y=`upper CI(95%)`+movement0),alpha=0.7, linetype=2) +
  geom_errorbar(data=m_edu2_icpp09, aes(x=x, ymin=lci+movement0,ymax=uci+movement0), alpha=0.8, width=2) + 
  geom_point(data=m_edu2_icpp09, aes(x=x, y=est+movement0), size=1.5, color="gray10", alpha=0.8) + 
  geom_text(data=m_edu2_icpp09, aes(x=x, y=-0.35+movement0, label=lb), family="serif") + 
  facet_grid(.~dv) + 
  scale_fill_manual(values=c("black","gray70")) + 
  scale_x_continuous(breaks=c(1975,1990,2000,2010)) +
  scale_y_continuous(labels = function(y) round(y - movement0,2)) +
  coord_cartesian(ylim=c(-0.04,0.95)) + 
  labs(x = "Year of Turning 19 (Typical Year of Entering University)",
       y = "Conditional Coefficient with 95% Confidence Interval",
       caption = paste0("Histogram indicates the frequency of cases (black = university educated, gray = not).")) +
  theme_classic(base_family="serif") + 
  theme(plot.title = element_text(hjust=0.5),
        panel.background = element_rect(color="black", size=1),
        strip.text = element_text(size=11, face="bold"),
        legend.position="none")
ggsave("mx_cplot_icpp09_alt.pdf", width=7, height=4)

#+
###################################################
## Analysis of 2013 Survey (Figure A3, Table A4) ##
###################################################

## Models with Categorical Interaction
mA_icpp13 <- lm_robust(immigincrease_alt~edu2*cohort+
                         male*cohort+
                         incomecat*cohort+
                         workstat*cohort+
                         married*cohort+
                         urban*cohort, 
                       data=icpp13, se_type="stata", clusters=as.factor(area))
mB_icpp13 <- lm_robust(foreignrights~edu2*cohort+
                         male*cohort+
                         incomecat*cohort+
                         workstat*cohort+
                         married*cohort+
                         urban*cohort, 
                       data=icpp13, se_type="stata", clusters=as.factor(area))
screenreg(list(mA_icpp13,mB_icpp13), 
          digits=3, include.ci=FALSE, single.row=TRUE,
          custom.model.names = c("Admit More Foreigners","Integrate Foreign Residents"),
          custom.coef.map = cfmap)
texreg(list(mA_icpp13,mB_icpp13), 
       digits=3, include.ci=FALSE, single.row=TRUE,
       custom.model.names = c("Admit More Foreigners","Integrate Foreign Residents"),
       custom.coef.map = cfmap, 
       custom.note = "%stars. Robust standard errors in parentheses.",
       booktabs = TRUE, dcolumn = TRUE, use.packages = FALSE, table=FALSE,
       file = "mtab_rc2_icpp13.tex")

mA_edu2_icpp13 <- 
  as.data.frame(  rbind(
    summary(lm_robust(immigincrease_alt~edu2*relevel(cohort,levels(cohort)[1])+
                        male*relevel(cohort,levels(cohort)[1])+
                        incomecat*relevel(cohort,levels(cohort)[1])+
                        workstat*relevel(cohort,levels(cohort)[1])+
                        married*relevel(cohort,levels(cohort)[1])+
                        urban*relevel(cohort,levels(cohort)[1]), 
                      data=icpp13, 
                      se_type="stata", clusters=as.factor(area)))$coefficients[2,c(1,4,5,6)],
    summary(lm_robust(immigincrease_alt~edu2*relevel(cohort,levels(cohort)[2])+
                        male*relevel(cohort,levels(cohort)[2])+
                        incomecat*relevel(cohort,levels(cohort)[2])+
                        workstat*relevel(cohort,levels(cohort)[2])+
                        married*relevel(cohort,levels(cohort)[2])+
                        urban*relevel(cohort,levels(cohort)[2]), 
                      data=icpp13, 
                      se_type="stata", clusters=as.factor(area)))$coefficients[2,c(1,4,5,6)],
    summary(lm_robust(immigincrease_alt~edu2*relevel(cohort,levels(cohort)[3])+
                        male*relevel(cohort,levels(cohort)[3])+
                        incomecat*relevel(cohort,levels(cohort)[3])+
                        workstat*relevel(cohort,levels(cohort)[3])+
                        married*relevel(cohort,levels(cohort)[3])+
                        urban*relevel(cohort,levels(cohort)[3]), 
                      data=icpp13, 
                      se_type="stata", clusters=as.factor(area)))$coefficients[2,c(1,4,5,6)],
    summary(lm_robust(immigincrease_alt~edu2*relevel(cohort,levels(cohort)[4])+
                        male*relevel(cohort,levels(cohort)[4])+
                        incomecat*relevel(cohort,levels(cohort)[4])+
                        workstat*relevel(cohort,levels(cohort)[4])+
                        married*relevel(cohort,levels(cohort)[4])+
                        urban*relevel(cohort,levels(cohort)[4]), 
                      data=icpp13, 
                      se_type="stata", clusters=as.factor(area)))$coefficients[2,c(1,4,5,6)]
  )  )
colnames(mA_edu2_icpp13) <- c("est","p","lci","uci")
mA_edu2_icpp13$dv <- "Admit More Foreigners"
# unique(sort(icpp13$univyr))
# median(1948:1974); median(1975:1989); median(1990:1999); median(2000:2008) 
mA_edu2_icpp13$x <- c(1961,1982,1994.5,2004)

mB_edu2_icpp13 <- 
  as.data.frame(  rbind(
    summary(lm_robust(foreignrights~edu2*relevel(cohort,levels(cohort)[1])+
                        male*relevel(cohort,levels(cohort)[1])+
                        incomecat*relevel(cohort,levels(cohort)[1])+
                        workstat*relevel(cohort,levels(cohort)[1])+
                        married*relevel(cohort,levels(cohort)[1])+
                        urban*relevel(cohort,levels(cohort)[1]), 
                      data=icpp13, 
                      se_type="stata", clusters=as.factor(area)))$coefficients[2,c(1,4,5,6)],
    summary(lm_robust(foreignrights~edu2*relevel(cohort,levels(cohort)[2])+
                        male*relevel(cohort,levels(cohort)[2])+
                        incomecat*relevel(cohort,levels(cohort)[2])+
                        workstat*relevel(cohort,levels(cohort)[2])+
                        married*relevel(cohort,levels(cohort)[2])+
                        urban*relevel(cohort,levels(cohort)[2]), 
                      data=icpp13, 
                      se_type="stata", clusters=as.factor(area)))$coefficients[2,c(1,4,5,6)],
    summary(lm_robust(foreignrights~edu2*relevel(cohort,levels(cohort)[3])+
                        male*relevel(cohort,levels(cohort)[3])+
                        incomecat*relevel(cohort,levels(cohort)[3])+
                        workstat*relevel(cohort,levels(cohort)[3])+
                        married*relevel(cohort,levels(cohort)[3])+
                        urban*relevel(cohort,levels(cohort)[3]), 
                      data=icpp13, 
                      se_type="stata", clusters=as.factor(area)))$coefficients[2,c(1,4,5,6)],
    summary(lm_robust(foreignrights~edu2*relevel(cohort,levels(cohort)[4])+
                        male*relevel(cohort,levels(cohort)[4])+
                        incomecat*relevel(cohort,levels(cohort)[4])+
                        workstat*relevel(cohort,levels(cohort)[4])+
                        married*relevel(cohort,levels(cohort)[4])+
                        urban*relevel(cohort,levels(cohort)[4]), 
                      data=icpp13, 
                      se_type="stata", clusters=as.factor(area)))$coefficients[2,c(1,4,5,6)]
  )  )
colnames(mB_edu2_icpp13) <- c("est","p","lci","uci")
mB_edu2_icpp13$dv <- "Integrate Foreign Residents"
mB_edu2_icpp13$x <- c(1961,1982,1994.5,2004)

m_edu2_icpp13 <- rbind(mA_edu2_icpp13,mB_edu2_icpp13)
m_edu2_icpp13$dv <- factor(m_edu2_icpp13$dv,levels=unique(m_edu2_icpp13$dv))
m_edu2_icpp13$lb <- factor(c("I","II","III","IV"), levels=c("I","II","III","IV"))
m_edu2_icpp13

## Interflex (without background variables)
set.seed(12345)
mxA_icpp13 <- interflex(Y = "immigincrease_alt", D = "edu2", X = "univyr",
                        Z = c("male", "incomecat", "workstat", "married", "urban"),
                        cl="area",
                        data = icpp13, na.rm=TRUE, estimator = "kernel",
                        full.moderate = TRUE, bw=5, #CI=FALSE,
                        nboots = 1000, parallel = TRUE, cores = 3, theme.bw=TRUE)
# mxA_icpp13$figure
set.seed(12345)
mxB_icpp13 <- interflex(Y = "foreignrights", D = "edu2", X = "univyr",
                        Z = c("male", "incomecat", "workstat", "married", "urban"),
                        cl="area",
                        data = icpp13, na.rm=TRUE, estimator = "kernel",
                        full.moderate = TRUE, bw=5, #CI=FALSE,
                        nboots = 1000, parallel = TRUE, cores = 3, theme.bw=TRUE)
# mxB_icpp13$figure

# save(mxA_icpp13, mxB_icpp13, file="mx_icpp13_alt.RData")
# load("mx_icpp13_alt.RData")

tmpA <- as.data.frame(mxA_icpp13$est.kernel$`1`)
tmpA$dv <- "Admit More Foreigners"
tmpB <- as.data.frame(mxB_icpp13$est.kernel$`1`)
tmpB$dv <- "Integrate Foreign Residents"
out <- rbind(tmpA,tmpB)
out$dv <- factor(out$dv,levels=unique(out$dv))

tmpA <- data.frame(x=mxA_icpp13$hist.out$mids,
                   y=c(mxA_icpp13$count.tr$`0`,mxA_icpp13$count.tr$`1`),
                   tr=rep(c("0","1"), each=length(mxA_icpp13$hist.out$mids)))
tmpA$dv <- "Admit More Foreigners"
tmpB <- data.frame(x=mxB_icpp13$hist.out$mids,
                   y=c(mxB_icpp13$count.tr$`0`,mxB_icpp13$count.tr$`1`),
                   tr=rep(c("0","1"), each=length(mxB_icpp13$hist.out$mids)))
tmpB$dv <- "Integrate Foreign Residents"
outh <- rbind(tmpA,tmpB)
outh$dv <- factor(outh$dv,levels=unique(outh$dv))
outh$tr <- factor(outh$tr, levels=c("1","0"))

#+ fig.width=7, fig.height=4
## Figure
movement0 <- 0.3
ggplot(out) + 
  geom_vline(aes(xintercept=1975), linetype=1, color="gray80") + 
  geom_vline(aes(xintercept=1990), linetype=1, color="gray80") + 
  geom_vline(aes(xintercept=2000), linetype=1, color="gray80") + 
  geom_vline(aes(xintercept=2010), linetype=1, color="gray80") + 
  geom_col(data=outh, aes(x=x,y=y*0.0015,fill=tr), width=0.8) + 
  geom_hline(aes(yintercept=0+movement0), linetype=2, color="gray50") + 
  geom_ribbon(aes(x=X,ymin=`lower CI(95%)`+movement0,
                  ymax=`upper CI(95%)`+movement0),alpha=0.3) +
  geom_line(aes(x=X,y=TE+movement0)) +
  # geom_line(aes(x=X, y=`lower CI(95%)`+movement0),alpha=0.7, linetype=2) +
  # geom_line(aes(x=X, y=`upper CI(95%)`+movement0),alpha=0.7, linetype=2) +
  geom_errorbar(data=m_edu2_icpp13, aes(x=x, ymin=lci+movement0,ymax=uci+movement0), alpha=0.8, width=2) + 
  geom_point(data=m_edu2_icpp13, aes(x=x, y=est+movement0), size=1.5, color="gray10", alpha=0.8) + 
  geom_text(data=m_edu2_icpp13, aes(x=x, y=-0.35+movement0, label=lb), family="serif") + 
  facet_grid(.~dv) + 
  scale_fill_manual(values=c("black","gray70")) + 
  scale_x_continuous(breaks=c(1975,1990,2000,2010)) +
  scale_y_continuous(labels = function(y) round(y - movement0,2)) +
  coord_cartesian(ylim=c(-0.04,0.95)) + 
  labs(x = "Year of Turning 19 (Typical Year of Entering University)",
       y = "Conditional Coefficient with 95% Confidence Interval",
       caption = paste0("Histogram indicates the frequency of cases (black = university educated, gray = not).")) +
  theme_classic(base_family="serif") + 
  theme(plot.title = element_text(hjust=0.5),
        panel.background = element_rect(color="black", size=1),
        strip.text = element_text(size=11, face="bold"),
        legend.position="none")
ggsave("mx_cplot_icpp13_alt.pdf", width=7, height=4)

#+
###################################################
## Analysis of 2022 Survey (Figure A2, Table A3) ##
###################################################

## Models with Categorical Interaction
mA_survey22 <- lm_robust(immigincrease_alt~edu2*cohort+
                           male*cohort+
                           incomecat*cohort+
                           workstat*cohort+
                           married*cohort+
                           urban*cohort, 
                         data=survey22, se_type="stata")
mB_survey22 <- lm_robust(foreignrights~edu2*cohort+
                           male*cohort+
                           incomecat*cohort+
                           workstat*cohort+
                           married*cohort+
                           urban*cohort, 
                         data=survey22, se_type="stata")
screenreg(list(mA_survey22,mB_survey22), 
          digits=3, include.ci=FALSE, single.row=TRUE,
          custom.model.names = c("Admit More Foreigners (Nation)","Integrate Foreign Residents"),
          custom.coef.map = cfmap)
texreg(list(mA_survey22,mB_survey22), 
       digits=3, include.ci=FALSE, single.row=TRUE,
       custom.model.names = c("Admit More Foreigners (Nation)","Integrate Foreign Residents"),
       custom.coef.map = cfmap, 
       custom.note = "%stars. Robust standard errors in parentheses.",
       booktabs = TRUE, dcolumn = TRUE, use.packages = FALSE, table=FALSE,
       file = "mtab_rc1_survey22.tex")

mA_edu2_survey22 <- 
  as.data.frame(  rbind(
    summary(lm_robust(immigincrease_alt~edu2*relevel(cohort,levels(cohort)[1])+
                        male*relevel(cohort,levels(cohort)[1])+
                        incomecat*relevel(cohort,levels(cohort)[1])+
                        workstat*relevel(cohort,levels(cohort)[1])+
                        married*relevel(cohort,levels(cohort)[1])+
                        urban*relevel(cohort,levels(cohort)[1]), 
                      data=survey22, se_type="stata"))$coefficients[2,c(1,4,5,6)],
    summary(lm_robust(immigincrease_alt~edu2*relevel(cohort,levels(cohort)[2])+
                        male*relevel(cohort,levels(cohort)[2])+
                        incomecat*relevel(cohort,levels(cohort)[2])+
                        workstat*relevel(cohort,levels(cohort)[2])+
                        married*relevel(cohort,levels(cohort)[2])+
                        urban*relevel(cohort,levels(cohort)[2]), 
                      data=survey22, se_type="stata"))$coefficients[2,c(1,4,5,6)],
    summary(lm_robust(immigincrease_alt~edu2*relevel(cohort,levels(cohort)[3])+
                        male*relevel(cohort,levels(cohort)[3])+
                        incomecat*relevel(cohort,levels(cohort)[3])+
                        workstat*relevel(cohort,levels(cohort)[3])+
                        married*relevel(cohort,levels(cohort)[3])+
                        urban*relevel(cohort,levels(cohort)[3]), 
                      data=survey22, se_type="stata"))$coefficients[2,c(1,4,5,6)],
    summary(lm_robust(immigincrease_alt~edu2*relevel(cohort,levels(cohort)[4])+
                        male*relevel(cohort,levels(cohort)[4])+
                        incomecat*relevel(cohort,levels(cohort)[4])+
                        workstat*relevel(cohort,levels(cohort)[4])+
                        married*relevel(cohort,levels(cohort)[4])+
                        urban*relevel(cohort,levels(cohort)[4]), 
                      data=survey22, se_type="stata"))$coefficients[2,c(1,4,5,6)]
  )  )
colnames(mA_edu2_survey22) <- c("est","p","lci","uci")
mA_edu2_survey22$dv <- "Admit More Foreigners (Nation)"
# unique(sort(survey22$univyr))
# median(1963:1989); median(1990:1999); median(2000:2009); median(2010:2022) 
mA_edu2_survey22$x <- c(1976,1994.5,2004.5,2016)

mB_edu2_survey22 <- 
  as.data.frame(  rbind(
    summary(lm_robust(foreignrights~edu2*relevel(cohort,levels(cohort)[1])+
                        male*relevel(cohort,levels(cohort)[1])+
                        incomecat*relevel(cohort,levels(cohort)[1])+
                        workstat*relevel(cohort,levels(cohort)[1])+
                        married*relevel(cohort,levels(cohort)[1])+
                        urban*relevel(cohort,levels(cohort)[1]), 
                      data=survey22, se_type="stata"))$coefficients[2,c(1,4,5,6)],
    summary(lm_robust(foreignrights~edu2*relevel(cohort,levels(cohort)[2])+
                        male*relevel(cohort,levels(cohort)[2])+
                        incomecat*relevel(cohort,levels(cohort)[2])+
                        workstat*relevel(cohort,levels(cohort)[2])+
                        married*relevel(cohort,levels(cohort)[2])+
                        urban*relevel(cohort,levels(cohort)[2]), 
                      data=survey22, se_type="stata"))$coefficients[2,c(1,4,5,6)],
    summary(lm_robust(foreignrights~edu2*relevel(cohort,levels(cohort)[3])+
                        male*relevel(cohort,levels(cohort)[3])+
                        incomecat*relevel(cohort,levels(cohort)[3])+
                        workstat*relevel(cohort,levels(cohort)[3])+
                        married*relevel(cohort,levels(cohort)[3])+
                        urban*relevel(cohort,levels(cohort)[3]), 
                      data=survey22, se_type="stata"))$coefficients[2,c(1,4,5,6)],
    summary(lm_robust(foreignrights~edu2*relevel(cohort,levels(cohort)[4])+
                        male*relevel(cohort,levels(cohort)[4])+
                        incomecat*relevel(cohort,levels(cohort)[4])+
                        workstat*relevel(cohort,levels(cohort)[4])+
                        married*relevel(cohort,levels(cohort)[4])+
                        urban*relevel(cohort,levels(cohort)[4]), 
                      data=survey22, se_type="stata"))$coefficients[2,c(1,4,5,6)]
  )  )
colnames(mB_edu2_survey22) <- c("est","p","lci","uci")
mB_edu2_survey22$dv <- "Integrate Foreign Residents"
mB_edu2_survey22$x <- c(1976,1994.5,2004.5,2016)

m_edu2_survey22 <- rbind(mA_edu2_survey22,mB_edu2_survey22)
m_edu2_survey22$dv <- factor(m_edu2_survey22$dv,levels=unique(m_edu2_survey22$dv))
m_edu2_survey22$lb <- factor(c("I & II","III","IV","V"), levels=c("I & II","III","IV","V"))
m_edu2_survey22

## Interflex (without background variables)
set.seed(12345)
mxA_survey22 <- interflex(Y = "immigincrease_alt", D = "edu2", X = "univyr",
                          Z = c("male", "incomecat", "workstat", "married", "urban"),
                          data = survey22, na.rm=TRUE, estimator = "kernel",
                          full.moderate = TRUE, bw=5, #CI=FALSE,
                          nboots = 1000, parallel = TRUE, cores = 3, theme.bw=TRUE)
# mxA_survey22$figure
set.seed(12345)
mxB_survey22 <- interflex(Y = "foreignrights", D = "edu2", X = "univyr",
                          Z = c("male", "incomecat", "workstat", "married", "urban"),
                          data = survey22, na.rm=TRUE, estimator = "kernel",
                          full.moderate = TRUE, bw=5, #CI=FALSE,
                          nboots = 1000, parallel = TRUE, cores = 3, theme.bw=TRUE)
# mxB_survey22$figure

# save(mxA_survey22, mxB_survey22, file="mx_survey22_alt.RData")
# load("mx_survey22_alt.RData")

tmpA <- as.data.frame(mxA_survey22$est.kernel$`1`)
tmpA$dv <- "Admit More Foreigners (Nation)"
tmpB <- as.data.frame(mxB_survey22$est.kernel$`1`)
tmpB$dv <- "Integrate Foreign Residents"
out <- rbind(tmpA,tmpB)
out$dv <- factor(out$dv,levels=unique(out$dv))

tmpA <- data.frame(x=mxA_survey22$hist.out$mids,
                   y=c(mxA_survey22$count.tr$`0`,mxA_survey22$count.tr$`1`),
                   tr=rep(c("0","1"), each=length(mxA_survey22$hist.out$mids)))
tmpA$dv <- "Admit More Foreigners (Nation)"
tmpB <- data.frame(x=mxB_survey22$hist.out$mids,
                   y=c(mxB_survey22$count.tr$`0`,mxB_survey22$count.tr$`1`),
                   tr=rep(c("0","1"), each=length(mxB_survey22$hist.out$mids)))
tmpB$dv <- "Integrate Foreign Residents"
outh <- rbind(tmpA,tmpB)
outh$dv <- factor(outh$dv,levels=unique(outh$dv))
outh$tr <- factor(outh$tr, levels=c("1","0"))

#+ fig.width=7, fig.height=4
## Figure
movement2 <- 0.1
ggplot(out) + 
  geom_vline(aes(xintercept=1975), linetype=2, color="gray80") + 
  geom_vline(aes(xintercept=1990), linetype=1, color="gray80") + 
  geom_vline(aes(xintercept=2000), linetype=1, color="gray80") + 
  geom_vline(aes(xintercept=2010), linetype=1, color="gray80") + 
  geom_col(data=outh, aes(x=x,y=y*0.0008,fill=tr), width=0.8) + 
  geom_hline(aes(yintercept=0+movement2), linetype=2, color="gray50") + 
  geom_ribbon(aes(x=X,ymin=`lower CI(95%)`+movement2,
                  ymax=`upper CI(95%)`+movement2),alpha=0.3) +
  geom_line(aes(x=X,y=TE+movement2)) +
  # geom_line(aes(x=X, y=`lower CI(95%)`+movement2),alpha=0.7, linetype=2) +
  # geom_line(aes(x=X, y=`upper CI(95%)`+movement2),alpha=0.7, linetype=2) +
  geom_errorbar(data=m_edu2_survey22, aes(x=x, ymin=lci+movement2,ymax=uci+movement2), alpha=0.8, width=2) + 
  geom_point(data=m_edu2_survey22, aes(x=x, y=est+movement2), size=1.5, color="gray10", alpha=0.8) + 
  geom_text(data=m_edu2_survey22, aes(x=x, y=-0.115+movement2, label=lb), family="serif") + 
  facet_grid(.~dv) + 
  scale_fill_manual(values=c("black","gray70")) + 
  scale_x_continuous(breaks=c(1975,1990,2000,2010)) +
  scale_y_continuous(labels = function(y) round(y - movement2,2)) +
  coord_cartesian(ylim=c(-0.02,0.55)) + 
  labs(x = "Year of Turning 19 (Typical Year of Entering University)",
       y = "Conditional Coefficient with 95% Confidence Interval",
       caption = paste0("Histogram indicates the frequency of cases (black = university educated, gray = not).")) +
  theme_classic(base_family="serif") + 
  theme(plot.title = element_text(hjust=0.5),
        panel.background = element_rect(color="black", size=1),
        strip.text = element_text(size=11, face="bold"),
        legend.position="none")
ggsave("mx_cplot_survey22_alt.pdf", width=7, height=4)

